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Abstract 

We consider the design of a pattern recognition that matches templates to images, 
both of which are spatially sampled and encoded as temporal sequences. The image 
is subject to a combination of various perturbations. These include ones that can be 
\ modeled as parameterized uncertainties such as image blur, luminance, translation, 

^ ' and rotation as well as unmodeled ones. Biological and neural systems require that 

these perturbations be processed through a minimal number of channels by simple 
adaptation mechanisms. We found that the most suitable mathematical framework to 
meet this requirement is that of weakly attracting sets. This framework provides us 
with a normative and unifying solution to the pattern recognition problem. We analyze 
the consequences of its explicit implementation in neural systems. Several properties 
inherent to the systems designed in accordance with our normative mathematical ar- 
■ gument coincide with known empirical facts. This is illustrated in mental rotation, 

t — ■ visual search and blur/intensity adaptation. We demonstrate how our results can be 

applied to a range of practical problems in template matching and pattern recognition. 
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1 Notational preliminaries 

X 

We define an image as a mapping So(x,y) from a class of locally bounded mappings S C 
Loo(f2 x . x fly), where Q x C IR, Q y C R, and L oo (0 a; x Q y ) is the space of all functions 
/ : Q x x Q y —> M such that ||/||oo = esssup{||/(x, y)\\, x G fl x , y G Q y } < oo. Symbols x, y 
denote variables on different spatial axes. The value of Sq(x, y) depends on the characteristic 
value of interest (e.g. brightness, contrast, color, etc.). We assume that an image can be 
described within the system as a set of a-priory specified templates, Si(x, y) G S, i G X C N, 
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where X is the set of indices of all templates associated with the image Sq(x, y) G S. Symbol 
X + is reserved for X + = X U 0. 

The solution of a system of differential equations x = f (x, t, 6, u(t)), u : M> — * M m , 
6 G M. d passing through point x at t — t will be denoted for t > t as x(t, x , t , 6, u), 
or simply as x(t) if it is clear from the context what the values of Xo, 6 are and how the 
function u(t) is defined. 

By L^JtojX], to > 0, T > t we denote the space of all functions f : M>o — * IR n such 
that ||f ||oo,[to,T] = esssup{||f(t)||,t G [t ,T]} < oo; ||f ||oo,[t ,T] stands for the L^[t ,X] norm 
of f(t). 

Let A be a set in M n and || • || be the usual Euclidean norm in M. n . By the symbol ||-||^ 
we denote the following induced norm: 

ll x IU = ^f{||x-q||} 
In case x is a scalar and A G M>o, notation ||x||a stands for the following 



\x U 



|x| — A, |x| > A 
0, \x\ < A 



2 Introduction 

Template matching is the oldest and most common method for detecting an object in an 
image. According to this method the image is searched for items that match a template. A 
template consists of one or more local arrays of values representing the object, e.g. intensity, 
color, or texture. Between these templates and certain domains of the image, a similarity 
value is calculated^ , and a domain is associated with a template once their similarity exceeds 
a given threshold. 

Despite the simple and straightforward character of this method, its implementation 
requires us to consider two fundamental problems. The first relates to what features should 
be compared between the image S (x,y) and the template Si(x,y), i G X. The second 
problem is how this comparison should be done. 

The normative answer to the question of what features should be compared invokes solv- 
ing the issue of optimal image representation, ensuring most effective utilization of available 
resources and, at the same time, minimal vulnerability to uncertainties. Principled solu- 
tions to this problem are well-known from the literature and can be characterized as spatial 
sampling. For example, when the resource is frequency bandwidth of a single measurement 



1 Traditionally a correlation measure is commonly used for this purpose 
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measurement mechanism, the optimality of spatially sampled representations is proven in 
Gabor's seminal work [9q. In classification problems, the advantage of spatially sampled 
image representations is demonstrated in [12] . In general, these representations are obtained 
naturally when balancing the system resources and uncertainties in the measured signal. A 



simple argument supporting this claim is provided in Appendix 1 

The simplest form of spatial sampling can be achieved by factorizing both the domain 
fl x x fl y of the image So and the templates Si, i G X into subsets: 

fl x x fl y = \^jfl x ,t x £ly,ti t G fit, fl x ,t ^= Qy,t ^= &y (1) 
t 

Factorization flTJ induces sequences {S it \, where S itt are the restrictions of mappings Si to 
the domains fl x> t x fl y j- These sequences constitute sampled representations of Si, i G X + 
(see Figure [[]). Notice that the sampled image and template representations {Si,t} are, 
strictly speaking, sequences of functions. In order to compare them, scalar values f(S itt ) 
are normally assigned to each Si t- Examples include various functional norms, correlation 
functions, spectral characterizations (average frequency or phase), or simply weighted sums 
of the values of S^t over the entire domain fl Xtt x fl Vj t- Formally, / could be defined as a 



2 Consider, for instance, a system which measures image Si(x,y) using a set of sensors {m-i, . . . ,m n }. 
Each sensor m, is capable of measuring signals within the given frequency band Aj at the location Xi in 
corresponding spatial dimension x. Then according to [9], sensor rrii can measure both the frequency content 
of a signal and its spatial location with minimal uncertainty only if the signal has a Gaussian envelope in x: 
Si(x,y) - e CT = r In other words, the signal should be practically spatially bounded. This implies that 

the image must be spatially sampled. 
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functional, which maps restrictions S^t into the field of real numbers: 



/ 




This formulation allows a simple representation of images and templates as sequences of 
scalar values {/(S^t)}, i G 2T + , t G Q t - We will therefore adopt this method here. 

The answer to the second question, that of how the comparison is done, involves finding 
the best possible and most simple way to utilize information provided by a given image 
representation, at the same time ensuring invariance to basic distortions. Despite the fact 
that considerable attention has been given to this problem, a principled and unified solution 
is not yet available. The primary goal of our current contribution is to present a unified 
framework to solve this problem for a class of systems of sufficiently broad theoretical and 
practical relevance. 

We consider the class of systems in which spatially sampled image representations are 
encoded as temporal sequences. In other words, parameter t in the notation f(Si t) is the 
time variable. This type of representation is frequently encountered in neuronal networks 
[T3] (see also references therein), and so such systems have a claim to biological plausibility. 
In addition, they enable a simple solution to a well-known dilemma. The dilemma is about 
whether comparison between templates and image domains should be made on a large, global, 
or on a small, local scale. The solution to this dilemma consists in temporal integration. Let, 
for instance, Q t = [fl,T], T 6 K>o- Then an example of a temporally-integral, yet spatially 
sampled, representation is: 



The temporal integral 4>i(t) contains both spatially local and global image characterizations. 
Whereas its time-derivative at t equals to f(Si tt ) and corresponds to spatially sampled, local 
representation S^t, the global representation 0;(T) equals to the integral, cumulative charac- 
terization of Si. An example illustrating these properties is provided in Figure [21 A further 
advantage of spatiotemporal representations (pi(t) is that they offer powerful mechanisms for 
comparison, processing and matching of (j>i(t), i G X. These mechanisms can generally be 
characterized in terms of dynamic oscillator networks which synchronize when their inputs 
are converging to the same function. 

Despite advantages such as optimality, simplicity and biological plausibility, there are 
theoretical issues which have prevented wide application of spatiotemporal representations 
to template matching. The most important issues, from the authors' viewpoint, are, first, 



= f f(Si,r)dr, te[0,T] 



i G T 



■+ 



(3) 
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Figure 2: Spatiotemporal image representation via spatial sampling and temporal integra- 
tion. Panel a contains original object, So; (x,y) marks a point on the image with respect to 
which the correlation is calculated; factorization of the domain Q x x Q y into ten nonintersect- 
ing subsets fl x x Q y = \Jj < L 1 £l x>tj xriy > t j - Panel b - templates Si, S2 and plots of (Si^Xx, y), 
ft (S2,tj)( x , y) ~~ the values of the normalized correlation between = Si(Q Xttj x fi^) and 
So(£l x ,tj x &y,tj)- Panel c - plots of the values of ([3D as a function of parameter t for templates 
Si (blue line) and S 2 (red line). 

how to achieve effective recognition in the presence of modeled disturbances, of which the 
most common ones are blur, luminance, and rotational and translational distortion. Second, 
how to take into account inevitable unmodeled perturbations. 

The first class of problems amounts to finding an identification/adaptation algorithm 
capable of reconstructing parameters of generally nonlinear perturbations. Currently avail- 
able approaches either are restricted to linear parametrization of disturbances, involve over- 
parametrization, or use domination feedback. However, linear parametrization is too re- 
stricted to be plausible, overparametrization is expensive in terms of the number of ad- 
justable units, and domination lacks adequate sensitivity. For these reasons these methods 
remain unsatisfactory. The second class of problems calls for procedures for recognizing an 
image from its perturbed temporal representation <fri(t). At this level the system is facing 
contradictory requirements of ensuring robust performance while being highly sensitive to 
minor changes in the stimulation. 

Both these problems are traditionally dealt with within the concept of Lyapunov-stable 
attractors. By allowing the system to converge on an attractor, it is possible to eliminate 
modeled and unmodeled distortions and thus, for instance, complete an incomplete pattern in 
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Figure 3: General scheme of a system for adaptive template matching using temporal codes. 
Level 1 contains the adaptive compartments. Its functional role is to ensure invariance to 
modeled uncertainties. Level 2 corresponds to the comparison compartments and consists 
of coupled nonlinear oscillators. Solid arrows represent the information flow in the system. 

the input [H [HI HSl UHl E2] . The advantage of these methods resides in the robustness inherent 
in uniform asymptotic Lyapunov stability. This advantage, however, comes at a cost: such 
systems are generally lacking in flexibility. Each stable attractor represents one pattern; 
but often an image contains more than one pattern. When the system is steered to one 
template, the other is lost from the representation. It would, therefore, be preferable to have 
a system that allows flexible switching between alternative patterns. Yet, the very notion 
of stable convergence to an attractor prevents switching and exploration across patterns. 
Furthermore, as we will show, for a class of the images with multiple representations and 
various symmetries globally stable solutions to the problem of invariant template matching 
may not even exist. 

We propose a unifying framework capable of combining robustness and flexibility. In 
contrast to common intuitions, which aim at achieving desired robustness by means of stable 
attractors, we advocate instability as an advantageous substitute. We demonstrate that a 
specific type of instability, the concept of weakly attracting sets, provides both the necessary 
invariance and flexibility. 

To illustrate these principles we designed a recognition system consisting of two major 
subcomponents (see Figure [3]). The first is an adaptive component in which information 
is processed by a class of spatiotemporal filters. These filters represent internal models of 
distortions. The models of most common distortions, including rotation, translation, and 
blur, are often nonlinearly parameterized. Until recently adequate compensatory mechanisms 
for nonlinear parameterized uncertainties where unexplored territory. In recent work [JT] we 
have shown that the problem of non- dominating adaptation could, in principle, be solved 
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within the concept of Milnor or weak, unstable attractors. Here we provide a solution to this 
problem that will enable systems to deal with specific nonlinearly parameterized models of 
distortions that are typical for a variety of optical and geometrical perturbations. 

The second component consists of a network of coupled nonlinear oscillators. These op- 
erate as coincidence detectors. Each oscillator in our system represents a Hindmarsh-Rose 
model neuron. These model neurons are generally believed to provide a good qualitative 
approximation to biological neuron behavior. At the same time they are computationally ef- 
fective in simulations [20]. For networks of these oscillators we prove, first of all, boundedness 
of the state of the perturbed solutions. In addition, we specify the parameter values which 
lead to emergence of globally stable invariant manifolds in the system state space. Although 
we do not provide explicit criteria for met a- stability in this class of networks, the conditions 
presented allow us to narrow substantially the domain of relevant parameter values in which 
this behavior is to be found. 

There is an interesting consequence to the unstable character of the compensation for 
modeled perturbations. When the system negotiates multiple classes of uncertainties simul- 
taneously (e.g. focal/contrast and intensity/luminance), different types of compensatory 
adjustments occur at different time scales. Adaptation at different time scales is a well- 
known phenomenon in biological visual systems, in particular when light/dark adaptation is 
combined with optical/neuronal blur [TH [23, [29], [33] ; experiments have shown that combined 
adaptation processes take place simultaneously but at different time-scales [3j [5j [3H [35] . Our 
analytical study suggests that this difference in time-scales emerges naturally as a sufficient 
condition for the proper operation of our system. 

This paper is organized as follows. In Section [3] we provide a formal description of the 
class of images and templates, and formally state the problems of our study. In Section [4] we 
provide the main results of our present contribution. In Section [5] we discuss the theoretical 
results, relate them to relevant observations in the empirical literature on visual perception 
and adaptation, and provide an application of our approach to a realistic pattern recognition 
problem: the detection of morphological changes in dendritic spines based on measurements 
obtained from multiphoton scanning microscope. 

3 Preliminaries and problem formulation 

We assume that the values Sq(x, y) of the original image are not available explicitly to the 
system; the system is able to measure only perturbed values of So(x,y). Perturbation is 
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defined as a mapping T: 

F[S Q , 0} : L^O, x Q y ) xl d -4 Looin* x fi„), 

where is the vector of parameters of the perturbation. The values of 6 are assumed to be 
unknown a-priori, whereas the mapping T is known. 

In systems for processing spatial information, mappings T often belong to a specific class 
that can be defined as follows: 

F[S , 0] = 0i ■ F[S , 9 2 ), 9 1 e R, 9 2 e R 

F[S , 6 2 ] : x Q y ) xR^ Loo(fi x x (4) 

O = (6 1 ,0 2 ) 

Parameter #1 G [0i, m i n , #i,max] C R in ((3]) models linear perturbations, for instance variations 
of overall brightness or intensity of the original image Sq. It can be interpreted also as an 
a-priori unknown gain in the measurement channel of a sensor. Mapping F(So,9 2 ) in Q, 
parameterized by 9 2 G [#2, m in, #2,max] C R, corresponds to typical nonlinear perturbations 
of image So- Table CD provides examples of these perturbations, their mathematical models 
and the physical meaning of parameter 9 2 . Throughout the paper we assume that mappings 
F[S , $2] are Lipschitz in 9 2 : 

3 D G R >0 : \f[S ,9' 2 ](x,y) - F[S ,9'£\(x,y)\ < D\9' 2 - 9' 2 % 

(5) 

V (x, y) G tt x x Q y , 9' 2 , 9' 2 ' G R 
Notice that, strictly speaking, several typical transformations such as translation, scaling, 
and rotation, are not always Lipschitz. This is because image So can, for instance, have sharp 
edges which corresponds to discontinuities in x, y. In practice, however, prior application 
of a blurring linear filter will render sharp edges in an image smooth, thus assuring that 
condition (jSJ) applied 

The image JF[So,0] is assumed to be spatially sampled according to factorization ([T|): 

HSo,0](x,y) = { ^ e x t G Q t (6) 

Because index t in ([6]) is assumed to be a time variable we let Q t — [0, 00). To each Ft [So, 0} 
a value /(Ft [So, 6]) G R is assigned. Formally this procedure can be defined by a functional 
which maps mappings .Ft [So, 0} into the real values: 

/ : LooiO, x n y ) -> R. (7) 

3 In biological vision discontinuity of So in x, y corresponds to images with abrupt local changes in bright- 
ness along spatial dimensions x, y. Although this is a rather common situation in nature, in visual systems 
actual images So rarely reach a sensor in their spatially discontinuous form. In fact, prior to reaching the 
sensory part, they are subject to linear filtering induced by optics. Therefore the images that reach the 
sensor are always smooth. Hence condition ([3]) will generally be satisfied. 
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Table 1: Examples of typical nonlinear perturbations of Sq. Parameter Ag in the right 
column is a positive constant 



Physical meaning 


Mathematical model 

ofns o ,0 2 ] 


Domain of 
physical relevance 


Translation (in x dimension) 
2 - shift 


nS o ,02] = S o (x + 02,y) 


-A fl < 2 < A e 


Scaling (in x dimension) 
2 - scaling factor 


nS o ,0 2 ] = S o (02-x,y) 


< 02 < A e 


Rotation 
around the origin 

82 - angle of rotation 


F[So, 02] = S (x r (x, y, 2 ), y r (x, y, 6 2 )) 

x r (x, y, 02) = cos(02)x - sm(0 2 )y 
y r (x, y,9 2 )= sin(0 2 )x + cos(0 2 )y 


< 02 < 2vr 


Image blur [5] 
(not normalized) 

02 - blur parameter 


1 J IjrdiUbblcLIl. 

h = exp-i«^ )2+( ^ )2) 
2) Out-of-focus: 

1 0, else 


< 02 < A e 



In the singular case, when Q xt x VL yt is a point (x t ,yt), the mapping !F t [So,d](x,y) and 
functional / will be defined as / = J-t[So, B](xt, yt) = J~[So, 0](x t , yt)- 

We concentrated our efforts on obtaining a principled solution to the problem of invariant 
template matching in systems with spatiotemporal processing of information. For this reason 
we prefer not to provide a specific description of functionals /. We do, however, restrict our 
consideration to linear and Lipschitz functionals, e.g. the functionals satisfying the following 
constraints: 

f(KT) = Kf(F), V«eE, \f(F)-f(r)\ < D2\\F-r\U D 2 eR >0 (8) 

Examples of functionals / satisfying conditions ([5} and their physical interpretations are 
provided in Table [21 
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Table 2: Examples of spatially-sampled representations of So 



Physical meaning 


Mathematical model of / 


Spectral power within 
the given frequency bands: 

L0 X € [u a ,Ub], UJy € [w c ,W(i] 




Weighted sum 

I tot limtflTirp ronA/'nlntin'n 

with exponential kernel) 


/ = /o :c xn !/ ^^o,0](x,y)e-l— °Hv-v°\dxdy 
(xo,yo) is the reference, "attention" point 


Scanning the image 
along a given trajectory 

(*(*)> = (£(<)> 7(f)) 


n*,t x n y , t = (£(*), 7(f)) 

/ = ^o,0](e(t), 7 (f)) 



Taking into account (HI), ([6]) and the fact that / is linear, the following equality holds 

f(F t [s ,o}) = M-^oAD, ?t[s ,e 2 ) = | F[SoM{*,y), M) e x (9) 

For the sake of compactness, in what follows we replace f(Ft[So, 62]) in the definition of 
f(J r t [S ,8]) in (jnj) with the following notation 

/(•Ft[Sb, #2]) = /o(t, f 2 ),/o:fl t xI^t (10) 

Notation fo(t, 62) in (TTOT) allows us to emphasize the dependence of / on unknown 62, time 
variable £, and original image Sq. Subscript "0" in fllOp indicates that /o(£, #2) corresponds 
to the sampled and perturbed So (equations (@J, ([7]), (JSj)), and argument 6*2 is the nonlinear 
parameter of the perturbation applied to the image. Adhering to this logic, we introduce 
the notation 

where subscript "i" indicates that fi{t,9o) corresponds to the perturbed and sampled tem- 
plate Si, and 6*2 is the nonlinear parameter of the perturbation applied to the template. 

Let us now specify the class of schemes realizing temporal integration of spatially sampled 
image representations. Explicit realization of temporal integration ([3]) is not optimal because 
it might lead to unbounded outputs for a wide class of relevant signals, for instance signals 
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that are constant or periodic with a nonzero average. The behavior of a temporal integrator 
([3]), however, can be fairly well approximated by a first-order linear filter. For the sampled 
image and template representations #2), these filters can be defined as follows: 

00 = — 0o + A; -0i./o OA) 

\ (11) 

01 = — <k + k ■ ejiit, 2 ), k,r e R >0 , i e I 

T 

In contrast to Q, for filters fill I) it is ensured that their state remains bounded for bounded 
inputs. In addition, on a first approximation, equations (JTTJ) present a simple model of neural 
sensors, collecting and encoding spatially-distributed information in the form of a function of 
timaj. With respect to the physical realizability of (llljl . in addition to requirements (J5j), (13) 
we shall only assume that spatially sampled representations 9ifi(t, 9 2 ), i € T + of Si ensure 
the existence of solutions for system (fTTT) . 

Consider the dynamics of variables (po(t) and (f>i(t), % 6 X defined by fill I) . We say that 
the i-th template matches the image iff for some given e G M>o the following condition holds 

Hmsup|0 o (t)-0i(t)| <e. (12) 

The problem, however, is that parameters 81, 82 in (TlTI) are unknown a-priori. While pertur- 
bations affect the image directly, they do not necessarily influence the templates. Rather to 
the contrary, for consistent recognition the templates are better kept isolated from external 
perturbations - at least within the time frame of pattern recognition, although they may, 
in principle, be affected by adaptive learning on a larger time scale. Having fixed, unmodi- 
fied templates in comparison with perturbed image representations implies that even in the 
cases when objects corresponding to the templates are present in the image, temporal image 
representation <fio(t) will likely be different from any of the templates, <pi(t). This will render 
the chances that condition ffT2l) is satisfied very small, so a template would almost never be 
detected in an image. 

We propose that the proper way for the system to meet requirement (fT2l) is to mimic the 
effect of disturbances in the template. In order to achieve this template matching system 
should be able to track the unknown values of parameters 6±, 62- Hence the original equations 



4 In principle, equation can be replaced with a more plausible model of temporal integration such 
as integrate-and-fire, Fitzhugh-Nagumo, or Hodgkin-Huxley model neurons. These extensions, however, are 
not immediately relevant for the purpose of our current study. Therefore for the sake of clarity we decided 
to keep the mathematical description of the system as simple as possible, keeping in mind the possibility of 
extension to a wider class of temporal integrators (JTTJ) . 
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for temporal integration (11 ip will be replaced with the following 

0o = --0o + k ■ 9 1 f (t, 9 2 ) 

\ , , ( 13 ) 

4>i = — 4>i + k ■ 9i tl fi(t,9i 2 ), k,reR >0 , iel 

T 

where 8^, 6^2 are the estimates of #1, 62 ■ The estimates 6n, 9 i2 must track instantaneous 
changes of 9\, 62- The information required for such an estimation should be kept at the 
minimal possible level. An acceptable solution would be a simple mechanism capable of 
tracking the perturbations from the measurements of the image alone. The formal statement 
of this problem is provided below: 

Problem 1 (Invariance) For a given image Sq, template Si, and their spatiotemporal rep- 
resentations satisfying (TJP, (EJ], and l[T3\). find estimates 

= 9i,l(t,T, K, 0o, 0i), 9 it 2 = 9i s2 (t,T,K,(/)o,(l>i) (14) 

as functions of time t, variables <po, (pi and parameters r, re such that for all possible values 
of parameters 61 G [6>i, m i n , 9 hmax ], 9 2 G [9 2 ,mirn "2, max I 

1 ) solutions of system [TB]) are bounded; 

2) in case /o = fi property ( TJj| ) is ensured, and 

3) the following holds for some 9[ G [6>i, mm , 0i, max ], 9' 2 G [6> 2 ,min, 02,max]-' 

limsup \9 ifl (t, r, re, o (t), (pi{t)) -9' iX \ < e g ,i, £e,i e 

. ' (15) 

limsup \9 ij2 (t, r, re, (p (*),&(*)) ~ ^,2! < ^0,2, ^,2 e M+ 

t— >oo 

Once the solution to Problem [1] is found, the next step will be to ensure that similarities 
(I12p are registered in the system. Following the spirit of neural systems and in agreement 
with the structure in Figure [31 we propose that detection of similarities is realized by a 
system of coupled oscillators. In particular, we require that states of oscillators % and 
converge as soon as the signals 0o(£), 0i(£) become sufficiently close. 

In the present article we restrict ourselves to the class of systems composed of linearly 
coupled Hindmarsh-Rose model neurons [17] . This choice is motivated by the fact that these 
oscillators can reproduce a broad class of behaviors observed in real neurons while being 
computationally efficient [2U]. A network of these neural oscillators can be mathematically 
described as follows: 

{j^ = -axf + bx\ + yt - Zi + I + Ui + 4>%{t), 
Vi = c-dx\- y h i G X + (16) 

^ =e{s(xi + x ) - z^, 
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Variables Xj, yi, Zi correspond to membrane potential, and aggregated fast and slow adap- 
tation currents, respectively. Coupling Ui in ( fl6i) is assumed to be linear and symmetric: 



u 



/ u \ 

"1 

\u n J 



( Xq \ 
\%n J 



1 



( —n 
1 

V 1 



-n 



1 \ 
1 



-n 



(17) 



/ 



and parameter 7 G R + . Our choice of the coupling function in (|T7j) is motivated by the 
following considerations. Fist, we wish to preserve the intrinsic dynamics of the neural 
oscillators when they synchronize, e.g. when x t = xj, yi = yj, Zi = Zj, i,j G {0, . . . , n}. For 
this reason it is desirable that the coupling vanishes when the synchronous state is reached. 
Second, we seek for a system in which synchronization between two arbitrary nodes, say 
the 2-th and the j-th nodes, is determined exclusively by the degree of (mis) matches in 
4>i(t), 4>j{t), independently of the activity of other units in the system. Third, the coupling 
should "pull" the system trajectories towards the synchronous state. Coupling function (fT7j) 
satisfies all these requirements. 

We set parameters of equations (Jl6]l to the following values: 

a = 1, b = 3, c=l, d = 5, , , 

s = A Xn = 1.6, e = 0.001, 1 } 



which correspond to the regime of chaotic bursting in each uncoupled element in ( 1161) 

The problem of detection of similarities in <j)o{t) and 4>i(t) can now be stated as follows. 



Problem 2 (Detection) Let system [W\) . JT7| ) be given and there exist i G X such that 
condition [W\l is satisfied. Determine the coupling parameter 7 as a function of system !jJJ)$ 
parameters such that 

1) solutions of the system are bounded for all bounded <f> i; i G X; 

2) states (xo(t), yo(t), zo(t)) and (xi(t),yi(t), Zi(t)) asymptotically converge to a vicinity 
of the synchronization manifold xq = Xi, yo = y*, Zq — z^. In particular, 

limsup |xo(t) — Xi(t)\ < 5(e) 

limsup \y (t) - yi(t)\ < 5(e) 

t^oo 

limsup |^o(^) — z i(t)\ ^ 3(. £ ): 



t— >oo 



where 5(-) is a non- decreasing function vanishing at zero. 

In the next section we present solutions to the problems of invariance and detection. 
We start from considerations of what would be the most adequate concept of analysis. Our 
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considerations will lead us to conclusion that for solving the problem of invariance, using the 
concept of Milnor attractors is advantageous over traditional concepts resting on the notion 
of Lyapunov stability. This implies that the sets to which the estimates 9^i, 6^2 converge 
should be weakly attracting rather than Lyapunov stable. We present a simple mechanism 
realizing this requirement for a wide class of models of disturbances. With respect to the 
second problem, the problem of detection, we provide sufficient conditions for asymptotic 
synchrony in system ([161) . 

4 Main Results 

Consider a system of temporal integrators, (TT3]) . in which the template subsystem (second 
equation in (fl3l) ) is designed to mimic the temporal code of an image using adjustment mech- 
anisms fll4p . Ideally, the template subsystem should have a single adjustment mechanism, 
which is structurally simple and yet capable of handling a broad class of perturbations. In 
addition it should require the least possible amount of a-priori information about images 
and templates. 

To search for a possible adaptation mechanism let us first explore the available theoretical 
concepts which can be used in its derivation. The problem of invariance, as stated in Problem 
[U can generally be understood as a specific optimization task. Particular solutions to such 
tasks as well as choice of the appropriate mathematical tools depend significantly on the 
following characteristics: uniqueness of the solutions, convexity with respect to parameters, 
and sensitivity to the input data (images and templates). Let us consider wether the invariant 
template matching problem meets these requirements. 

Uniqueness. Solutions to the problem of invariant template matching are generally not 
unique. The image may contain multiple instances of the template. Even if there is only 
a single unique object the template may fit it in multiple ways, for instance because it has 
rotational symmetry. Both cases are illustrated in Figure HI A similar argument applies to 
translational invariance in the images with multiple instances of the template (right picture 
in Figure HJ). 

Non-linearity and non- convexity. The problem of invariant template matching is gener- 
ally nonlinear and nonconvex in 9\, 6 2 . The nonlinearity is already evident from Table [TJ 
To illustrate the potential non-convexity consider, for instance, the following function 

9 1 f i (t,6 2 ) = 6 1 [ e -l*-*o|-|y-wl j j e ~i- 2 ({-i) 2 +(y-i? ^(£,7)^7) dxdy (19) 

Jn x . t xn yt t \Jn x xn y J 
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template Image 1 Image 2 




Figure 4: Example of a template and images which lead to non-unique solutions in the 
problem of invariant template matching. Image 1 is a rotated version of the template. 
Because the template has rotational symmetry, the angles 9 2 = 9 2 ± § n, n = 0, 1, . . . at 
which the template and the image match to each other are not unique. Image 2 contains 
two multiple instances of the template, which also leads to non-uniqueness. 



which is a composition of the Gaussian blur model (the forth row in Table [T]) with spatial 
sampling and subsequent exponential weighting (the second row in Table [2]) . In the literature 
on adaptation two versions of the convexity requirement are available. The first version 
applies to the case where the difference 6ifi(t, 9 2 ) — 0i,ifi{t, 9i j2 ) is not accessible for explicit 
measurement, and the variables <fio(t), 4>i(t) should be used instead. In this case the convexity 
condition will have the following form [7J: 



1 ~~ ^wT^ — ^1,1/1(^)^1,2) 



difi(t,9 2 ) 







i,2j 



de. 



■&i,lfi(t, #1,2) 



> 



(20) 



9i,lfi(t, 0i,2) 



Term ej(0o, <f>i) in ([20]) is usually the difference ei(4>o,<f>i) = <fio ~ & an d has the meaning 
of error. For the same pairs of points 6 I 1 ,6' 2 and 9n, 9 i 2 condition fl20l) may hold of fail 
depending on the sign of ej(0o(i) , 4>i{t)) at the particular time instance t. Hence it is not 
always satisfied, not even for convex ^,1/1(^,^,2)- 

The second version of the convexity requirement applies when the difference 9ifi(t, 9 2 ) — 
9i,ifi{t,9ip) can be measured explicitly. In this case the condition is formulated as definite- 
ness of the Hessian of 9ifi(t,9 2 )- It can easily be verified, however, that this depends, for 
instance, on the values of Si (£,7) in ffl9l) . Hence both versions of the convexity conditions 
generally fail in invariant template matching. 

Critical dependence on stimulation. An important feature of of invariant template match- 
ing problem is that its solutions critically depend on particular images and templates. Pres- 
ence of rotational symmetries in the templates affect the number of solutions. Hence objects 
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with different number of symmetries will be characterized by sets of solutions with different 
cardinality. 

We conclude that the problem of invariant template matching generally assumes multiple 
alternative solutions, nonlinearity and non-convexity with respect to parameters, and the 
structure of its solutions depends critically on a-priori unknown stimulation. What would 
be a suitable way to approach this class of problems in a principled manner? 

Traditionally, processes of matching and recognition are associated with convergence of 
the system's state to an attracting set. In our case the system's state is defined by vector x: 

x — (00, 01) ■ ■ ■ j 0ij • • • j #1, 1, ^2,1, ■ ■ ■ #i,2> • • • ) 

The attracting set, A, is normally understood as a set satisfying the following definition [12]: 



Definition 1 A set A is an attracting set iff it is 

i) closed, invariant, and 

ii) for some neighborhood V of A and for all x G V the following conditions hold: 

x(t,xo) e V V t > 0; (21) 
lim ||x(t,Xo)|Lt = (22) 

Traditional techniques for proving attractivity employ the concept of Lyapunov asymptotic 
stability^. Although the notion of set attractivity is wider, the method of Lyapunov func- 
tions is constructive and, in addition, Lyapunov asymptotic stability implies the desired 
attractiviy. For these reasons it is highly practical, and the tandem of set attractivity in 
Definition [1] and Lyapunov stability has been used extensively in recognition systems, in- 
cluding Hopfield networks, recurrent neural nets, etc. 

The problem of invariant template matching, however, challenges the universality of these 
concepts. First, because of inherent non-uniqueness of the solutions, there are multiple 
invariant sets in the system's state space. Hence, global Lyapunov asymptotic stability 
cannot be ensured. Second, when each solution is treated as a locally stable invariant set, 
it is essentially important to know the domain of its attractivity. This domain, however, 
depends on properties of function 6ifo(t, 62) in ( fl"3l) which vary with stimulation. Third, no 
method exists for solving Problem [1] for general nonlinearly parameterized #1/0 (^$2) that 
assures Lyapunov stability of the system. 



5 We recall that the set A is (globally) Lyapunov asymptotically stable iff for all e > there exists 
5(xq,s) > such that ||xo|L < 5(x ,e) => ||x(f, x )||^ < e for all t > 0, and lim^oo ||x(t, xq)!^ = 
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Figure 5: Standard stable attractors, panel a, vs weak attractors, panel b. Domains of stable 
attractors are neighborhoods containing A±, A<i- Estimates of sizes of these domains depend 
on particular images So,i> So^, So,3- These estimates are depicted as closed curves around 
Ai, Ai- Once the state converges to either of the attractors it stays there unless, probably, 
when the image changes. In contrast to this, domains of attraction for Milnor attracting sets 
are not neighborhoods. Hence, even a slightest perturbation in the image induces a finite 
probability of escape from the attractor. Hence multiple alternative representations of the 
image could eventually be restored. 

In order to solve the problem of invariant template matching we therefore propose to 
replace the standard notion of attracting set with a less restrictive concept. In particular we 
use the concept of weak or Milnor attracting sets [28J: 

Definition 2 A set A is weakly attracting, or Milnor attracting set iff 

i) it is closed, invariant and 

ii) for some set V (not necessarily a neighborhood of A) with strictly positive measure 
and for all x G V limiting relation |HP holds 

The main difference between the notions of a weak attracting set, Definition [21 and 
the standard one, Definition [1], is that the domain of attraction is not required to be a 
neighborhood of A. On the one hand, this allows to us use mathematical tools beyond the 
concept of Lyapunov stability in order to avoid problems with nonlinear parametrization 
and critical dependance on stimulation. On the other hand, it offers a natural mechanism 
for systems to explore multiple image representations. This is illustrated in Figure [51 

In the next paragraph we present technical details of how Problem [T] could be solved 
within the framework of Milnor attractors. 
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4.1 Invariant template matching by Milnor attractors 

Consider system ( 1131) : 

<j>o = —<h + k-6 1 Mt,e 2 ) 

T 

fa = — fa + k ■ 9nfi(t, 8 i2 ), k,r e M>o, i g I 

T 

and assume that the i-th template is present in the image. This implies that both the image 
and the template will have, at least locally in space, sufficiently similar spatiotemporal 
representations. Formally this can be stated as follows: 

3Agr >0 : \eif (t,e 2 )-9 1 f i (t t e 2 )\ < a, v^a, t >o (23) 

Hence without loss of generality we can replace equations (fT3"|) with the following 

0o = — fa + k-6 1 f i (t,9 2 ) + e(t) 

fa = — & + fc-0i,i/i(*,0i,2), k,reR >0 , iEl 

T 

where e{t) G Loo[0, oo], ^(t)!^ < A is a bounded disturbance. Solving Problemd], therefore, 
amounts to finding adjustment mechanisms (1141) such that trajectories 4>o(t), 4>i(t) in (1241) 
converge and limiting relations ()15p hold. 

The main idea of our proposed solution to this problem can informally be summarized 
as follows. First, we introduce an auxiliary system 

A = g(A, 0o, fat), A G M A , g:K A xlxlxl> ^ R A (25) 

and define 8i t i, 8^ 2 as functions of A, 0o, and fa: 

Oi,i = #o(A, r, k, 0o, fa), &i,2 = 9 it2 (X, r, k, fa, fa). (26) 

Second, we show that for some e G M>o, and Q\ C M A the following set 

Q* = {fa, fa g M,A G K A | |0 o (t) -fa(t)\ <e, XeQ x C M A } 

is forward- invariant in the extended system ( 1241) . ( 1251) and ( 1261) . Third, we restrict our 
attention to systems which have a subset Q in their state space such that trajectories starting 
in Q converge to Q*. Finally, we guarantee that the state will eventually visit domain Q thus 
ensuring that (fl2l holds. 



18 



We have found that choosing extension (1251) in the class of simple third-order bilinear 

systems 

7i 

Ai = — • (0o - 4>i) 

T 

A 2 = 72 • A 3 • ||0o - 0i|U 7i,72 GK>o (27) 

A 3 = -72 • A 2 ■ ||0o -(pi\\e, \Ai(to) + A|(* ) = 1 
ensures solution to Problem [U Specific technical details and conditions are provided in 
Theorem [T] 

Theorem 1 Let system (24\ ), l^2l\ ) be given, and function fi(t,8 2 ) be separated from zero and 
bounded. In other words, there exist constants D 3 , G M >0 such that for all t > 0, 9 2 G 
[^2,min, ^2,max] following condition holds: 

D 3 < fi(t,9 2 ) < D 4 (28) 



Taen there exist positive 71, 72, and e (see Table\3\for the particular values): 

D 

D 3 



72 < 7i, £ > ( 1 + -± ) (29) 



sitcn i/iai adaptation mechanisms 

= + Ai 

— ^2,min + (A 2 (t) 



# 2 ,max — # 2 ,min 

(30) 



2 

deliver a solution to Problem^ In particular, for all Q\ G [0i, m in> 0i,max]> #2 G [02,minj $2,r 
the following properties are guaranteed: 



limsup |0 o (t) - 0i(t)| < e\ 3 9' 2 E [9 2Ma , # 2 , ma x] ■ lim ij2 (t) = # 2 



t— »oo 



t— >oo 



2- 



where the value of e, depending on the choice of parameters 72, 71, can &e made arbitrarily 
close to tA (1 + D4/D3). 



Proof of the theorem is provided in Appendix 2 



Let us comment on the conclusions and conditions of Theorem [U First of all, the theorem 
shows that each z-th subsystem ensuring invariance of spatiotemporal image representation 
to the given modelled perturbations can be composed of no more than four differential 
equations: 

Temporal integration : 0j = 0, + k ■ §nfi(t, 0j 2 ) (31a) 

T 

7i 

Fast adaptation dynamics : Ai = — • (0o — 0j) (31b) 

T 

Slow adaptation dynamics : If 2 ^ 2 ,u 1 11 (31c) 

I A3 = -7 2 • A 2 • ||0o - (pi\\e 
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Table 3: Parameters of the compensatory mechanisms (1301) 



Parameter 


Values 


7i 


72 


£ 


£>T l A V 1+ W + 7 1 [ W ^ 
Mi = A + &0i !max Z>D 2 |0 2 , ma 


M -^A #2,max ~ #2,minl \ 

v aJ 2 j; 

x ^2, min| 


72 


7 2 <(^) 2 [^i, m ^2 (l + ^)( 


#2,max — #2,min \ 1 

2 



Notice that the time scales of temporal integration (13 lap , adaptation to linearly pa- 
rameterized uncertainties, ( I31bl) . and adaptation to nonlinearly parameterized uncertainties, 
(13 left , are different. Because of this difference in the time scales, subsystem (131bp is referred 
to as slow adaptation dynamics and subsystem (13 left as fast adaptation dynamics. The dif- 
ference between the time scales determines the degree of invariance and precision in the 
resulting system. For instance, as follows from Table El ratio 72/71 affects the value of e. 
This value defines the acceptable level of mismatches between an image and a template. In 
other words, it regulates the sensitivity of the system. The smaller the ratio 72/71, the higher 
the sensitivity. Ratio 72/(1/1") (see proof for details) affects the conditions for convergence. 

Slow adaptation dynamics, (131 cj) . can be interpreted as a searching, or wandering dynam- 
ics in the interval [#2,min, #2,max]- Its functional purpose is to explore the interval [# 2 ,min, #2,max] 
for possible values of ^.2 when models of perturbation are inherently nonlinear and no other 
choice except of explorative search is available. Solutions of the searching dynamics in (I31cl) 
are harmonic signals with time-varying frequency 72 1| — <M0lle- The larger the error, 
the higher the frequency of oscillation. When 7 2 ||0o(£) ~ 0i(^)l|e is constant, for instance 
equals to unit, equations (I31cl) reduce to 

^2 = A3 

(32) 

A3 — — A2 
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In general, every subsystem 

A 2 = 0*(A 2 , A 3 ,*), 

(33) 

A3 = 93(^2, A 3 ,t), 9-2,5-3 eC° 

generating dense trajectories A 2 (t) in [02, m in, 02,max] for some initial conditions A 2 (to), ^3(^0) 
and, at the same time, ensuring boundedness of A 2 (t), A 3 (t) for all t G M>o could be a 
replacement for (1321) in ( 13 left (see also [38]). Conclusions of the theorem in this case will 
remain the same except, probably, with respect to the choice of the particular values of 
7x, 72, £ in Table [3j Our present choice of subsystem (|32|) in (13 left as a prototype for the 
searching trajectory was motivated primarily by its simplicity in realization and linearity in 
state. 

The fast adaptation dynamics, (131bj) . corresponds to exponentially stable mechanisms. 
This can easily be verified by differentiating the difference #j,i(£) — Qi with respect to time 



(see also (H6l) in Appendix 2). The function of the fast adaptation subsystem is to track 
instantaneous changes in 61 exponentially fast in such a way that the difference Oi,i(t) — Q\ 
is determined mostly by mismatches #i, 2 (t) — # 2 . 

The problem of template matching is solved through the interplay of searching dynamics 
@i,2(t) — 6 2 (see Figure [6]) and the contracting dynamics expressed by (fio(t) — (fii(t). We use 
the results from [38] to prove the emergence of weakly (Milnor) attracting sets in the system 
state space. 

In principle, linearity of the uncertainty models in 6\ is not necessary to guarantee 
exponential stability of — &i- As has been shown in [3D], exponential stability of 

9i,i{t) — Q\ can be ensured by the same function #j,i(t) as in (1261) if we replace Oifi(t, 62) with 
fi{t, 61,62) '■ IR>o x 1 x 1 -* 1. Nonlinearities fi(t, 61,62), however, should be monotone in 
6\. In this case condition (1281) is to be replaced with the following 

n <f MMvlA) - fi(t, 61, 62) „ r „ „ , ( s 

^3 < ^ a - 4 ' V 2 G r2,min, ^maxj (34) 

6i,l — 6\ 

The general line of the proof will remain unaffected by this extension. 

The proposed compensatory mechanisms (|24|) . (127|) (|30|) are nearly optimal in terms of 
the dimension of the state of the whole system. Indeed, in order to track uncertain and 
independent 6\, 62 two extra variables are to be introduced. This implies that the minimal 
dimension of state of a system which solves Problem [1] equals three. Our four- dimensional 
system is therefore close to the optimal configuration. Furthermore, as follows from the 
proof of the theorem, the dimension of the slow subsystem could be reduced to one. Thus, 
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Figure 6: The interplay between temporal integration, (13 lap , and fast and slow adaptation 
(131bp . (I31cl) in the proposed solution to the problem of invariant template matching. Panel 
a. Contracting dynamics corresponds to the processes of temporal integration of a template 
and adaptation to linearly parameterized uncertainties. Searching dynamics is due to the 
adaptation to nonlinearly parameterized uncertainties. Panel b. Diagram of the phase 
portrait of system (I31al) . (131bl) . (I31cl) . Interaction between searching and contracting sub- 
systems forms a weakly attracting invariant set A. Its basin of attraction is not necessarily 
a neighborhood of A. This means that some trajectories starting in a small vicinity of A 
might eventually leave its neighborhood (dashed trajectory), while trajectories starting far 
away from A might enter such neighborhoods and remain there (solid line). 

in principle, a minimal realization could be achieved. In this case, however, boundedness of 
the state for every initial condition is no longer guaranteed. 

Theorem [T] establishes conditions for convergence of the trajectories of our prototype 
system ( 1241) . (|27|) (1301) to an invariant set in the system state space. In particular, when 
matching condition (1231) is met, it assures that temporal representation <pi(t) of the template 
tracks temporal representation ^o(^) of the image. In the next subsection we discuss how the 
similarity between these temporal representations can be detected by a system of coupled 
spiking oscillators. 

4.2 Conditions for synchronization of coincidence detectors 

Consider coincidence detectors (|T6l) . (flT|) . (|T8l) modeled by a system of coupled Hindmarsh- 
Rose oscillators. The goal of this section is to provide a constructive solution to Problem 
[2J First, we seek for conditions ensuring global exponential stability of the synchronization 
manifold of <po(t) = 0»(t) when <pi{t) are identical for each i. We do this by showing that solu- 
tions of the system are globally bounded, and for each pair of indexes i, j G {0, . . . , n} there 
exists a differentiable positive definite function V(xi,yi, Zi,Xj,yj, Zj), dV/dxi = —dV/xj 
such that V grows towards infinity with distance from the synchronization manifold and for 
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all bounded continuous <pi(t) = <pj(t) the following holds: 

V < -aV, a G R >0 . (35) 

When 4>i{t) ^ 4>j(t) equation (1351) implies that 

dV 

V < -aV + - </>&)). (36) 

Then using (I36I) and comparison lemma [23] we show that convergence of <pi(t) to 4>j(t) at t — * 
oo implies convergence of variables Xi(t), yi(t), Zi(t), Xj(t), Vj(t), Zj(i) to the synchronization 
manifold. The formal statement of this result is provided in Theorem [2] 



Theorem 2 Let system ^B\} be given, function u be defined as in [Ffy and functions 4>i(t), 
i £ {0, . . . ,n} be bounded. Then 

1) solutions of the system are bounded for all 7 G M + ; 

2) if, in addition, the following condition is satisfied 

7> 7 \, it + A, (37) 



(n+1) -a \ 2 
then for all i, j G {0, . . . , n} condition 



limsup \4>i(t) — 4>j(t)\ < e 

t^oo 



implies that 



limsup \xi(t) — Xj(t)\ < 5(e), 

t— >oo 

lim sup \yi(t)- yj (t)\<6(e), (3 j 



limsup |zj(t) — ^(t)| < 5(e). 
where 5 : M + — > M + a monotone and vanishing at zero function. 

Theorem [2] specifies the boundaries for stable synchrony in the system of coupled neural 
oscillators f|T6|) as a function of the coupling strength, 7, and parameters a, b, and d of 
a single oscillator. The last three parameters represent properties of the membrane and 
combined with xq, e, s and / completely characterize the dynamics of a single model neuron 
[TT] . ranging from single spiking to periodic or chaotic bursts. 

The distinctive feature of Theorem [2] is that it is suitable for analysis of systems with 
external time- dependent perturbations <f>i(t). This property is essential for the comparison 
task, where the oscillators are fed with time- varying inputs and the degree of their mutual 
synchrony is the measure of similarity between the inputs. 
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While the theorem provides us with conditions for stable synchrony, it allows us to 
estimate the domain of values of the coupling parameter 7 corresponding to potential inter- 
mittent, itinerant |22[ [23] or meta-stable regimes. In particular, as follows from Theorem [21 
a necessary condition for unstable synchronization in system (1161) is 



Notice that conditions (1391 . f[371) do not depend on the "bifurcation" parameter / which 
usually determines the type of bursting in the single oscillator. They also do not depend 
on the differences in the time scales defined by parameter e between the fast x, y, and 
slow, z, variables. Hence these conditions apply in a wide range of system behavior on the 
synchronization manifold. This advantage also has a downside, because conditions f[3"9"l) . (I3T1) 
are too conservative. This, however, seems to be a reasonable price for invariance of criteria 
f l39|) . fl371) with respect to the full range of dynamical behavior of a generally nonlinear 
system. 

5 Discussion 

We provided a principled solution to the issue of invariance in the problem of template 
matching. The pattern recognition problem can be solved using a network of nonlinear 
oscillators which synchronize when mismatches in the temporal representations of image and 
templates converge to zero. Although the solution to the latter problem is not normative we 
tried to keep the number of relevant parameters at minimum. In particular the dimension 
of the state of a single adaptation compartment is three which is minimal for generation of 
spikes ranging from periodic to chaotic bursts. Moreover, conditions (|39|) . (1371) allow us to 
choose coupling strength 7 as single control parameter for regulating stability/instability of 
the synchronous activity in the network. 

In this section we provide further extensions of the basic results of Theorems (H El discuss 
possible links between the normative part of our theory and known results in vision, and 
provide simple illustrations how particular systems for invariant template matching can be 
constructed using these results. 

5.1 Extension to the frequency-encoding schemes 

For the sake of notational simplicity we restricted our attention to temporal representations 
©, (Q of spatially sampled images. These encoding schemes can be interpreted as scanning 




(39) 
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Figure 7: Panel a. Temporal representation of a spatially distributed stimulus using fre- 
quency encoding. A stimulus (upper row) S(x) is spatially sampled by partitioning its 
domain into the union of intervals Q{. For each Qi an integral /, = f(J r i) = J n _ S(x)dx is 
calculated and a frequency uji is assigned. The resulting temporal representation (lower row) 
is expressed as the sum of two amplitude-modulated harmonic signals of frequencies, 0U3, 
Panel b. Temporal representation of a two-dimensional pattern. The pattern consists 
of black filled circles. The image domain is partitioned into a collection of horizontal and 
vertical strips. Dark domains correspond to higher frequencies. 

of an image over time. Yet, the results of Section H] apply to a broader class of encoding 
schemes. One example is frequency-coding used in many neural systems. Let us consider 
factorization ([6]) where in the notation J r t [So,0](x,y) symbol t is replaced with v. In order 
to extend the initial encoding scheme to the domain of frequency/spike rate encoding we 
introduce an additional linear functional f w as follows: 

Ut, F V [S Q , 0}) = J2 ■ t) ■ f{Fu[S , 0}), (40) 

V 

where h : R — > R is a bounded periodic function, and u u are distinct real numbers indexed 
by v . Function h[io v ■ t) in (140]) serves as a basis or carrier function generating periodic 
impulses of various frequencies u u . Thus each z/-th spatial sample of the image is assigned 
a particular frequency, and the amplitude of the oscillation is specified by /(^v[»So, 0\). 
Temporal representation of a one-dimensional stimulus according to encoding scheme (|40l) 
is illustrated in Figure [TJ panel a. 

This encoding scheme is plausible to biological vision, when frequencies u u are ordered 
according to relative position of domains £l x ,u, Qy,v to the center of the image. This cor- 
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responds, in particular, to the receptive fields in cat retinal ganglion cells [5|. Because the 
functional f w is linear in f(J-„[So,0]) and function h[io v ■ t) is bounded for all t, condition 
(JHJ) will be satisfied for f u . Hence the conclusions of Theorem [1] apply to these systems. 

5.2 Multiple representations of uncertainties 

Another property of system (|24p . (1271) . fl30l) . in addition to its ability to accommodate relevant 
encoding schemes such as frequency/rate and sequential/random scanning, is that each single 
value of 6 2 G (#2,mm, #2,max) induces at least two distinct attracting sets in the extended space. 
Indeed 

Xl(t) + A§(t) = const = 1 

along the trajectories of ( 1241) . ( |27|) . ( l30i) (see also the proof of Theorem CD). Hence for almost 
every value of A2 (except when A2 = ±1) in the definition of 02(t) in (J30l) there will always 
be two distinct values of A3: 



These give rise to distinct invariant sets in the system state space for the single value of 
62. The presence of two complementary encodings for the same figure is a plausible as- 
sumption that has been used in the perceptual organization literature to explain a range of 
phenomena, including perceptual ambiguity, modal and amodal completion, etc. See [15] . 
[T5] for a review. A consequence of the presence of multiple attractors corresponding to the 
single value of perturbation is that the time for convergence (the decision time) may change 
abruptly with small variations of initial conditions. The latter property is well documented 
in human subjects [UJ. Furthermore, the presence of two attractors with different basins for 
a single value of perturbation will lead to asymmetric distributions of decision times, which 
is typically observed in human and animal reaction time data [37] . 

5.3 Multiple time scales for different modalities in vision 

An important property of the proposed solution to the problem of invariance is that the time 
scales of adaptation to linearly and nonlinearly parameterized uncertainties are substantially 
different. This difference in time scales emerged naturally in the course of our mathematical 
argument as a consequence of splitting the system dynamics into a slow searching subsystem 
and a fast asymptotically stable one. This allowed us to prove emergence of unstable yet 
attracting invariant sets thus ensuring existence of a solution to the problem of invariant 
template matching.In particular, Theorem [T] requires that the time constants of adaptation 
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to linearly parameterized uncertainties (for instance, unknown intensity of the image) are to 
be substantially smaller than the time constants of adaptation to nonlinearly parameterized 
uncertainties (image blur, rotation, scaling etc.). Furthermore, as follows from Table [3j the 
larger the difference in the time scales the higher the possible precision and the smaller the 
errors. 

Even though the difference in time scales was motivated purely by theoretical consid- 
erations, there is strong evidence that biological systems adapt at different time scales to 
uncertainties from different modalities. For example, the time scale of light adaptation is 
within tens of milliseconds [15] while adaptation to "higher-order" modalities like rotation 
and image blur extends from hundreds of milliseconds to minutes |44j . In motor learning 
the evidence of presence of slow and fast adaptation at the time scale minutes is reported in 
[36] . These findings, therefore, motivate our belief that system ( |24|) . (I2T|) . (!30|) could serve 
as a simple, yet qualitatively realistic, model for adaptation mechanisms in vision, motor 
behavior, and decision making. 

5.4 Rot at ion- invariant matching and mental rotation experiments 

Let us illustrate how the results of Theorems [U [2] can be applied to template matching when 
an object is rotated over an unknown angle and its brightness is uncertain a-priori. In order 
to obtain a temporal representation of the image we use the frequency-encoding scheme (1401) 
as is illustrated in Figure [7J panel b. In particular we use the following transformation 



where 82 is the rotation angle of image Si(x,y) around its central point, 81 is the image 
brightness, function h{oj u ■ t) = sin 2 ^ • t), and 



is simply an integral of the rotated image Si by the angle 82 over the strip VL XyU x VL y ^ v . 

According to (T2"4"|) . f[2"Tl) . f[3"Uj) . fTTB]) the recognition system (see Figure E] for its general 
structure) can be described by the system of differential equations provided in Table HI 
Details of their implementation, specific values of the parameters, initial conditions, and the 
source files of a working MATLAB Simulink model can be found in [39J . 

We tested system performance for a variety of input images, in particular the class of 
Garner patterns [10] (see Figure the first row) . These patterns are widely used in the 
experiments with humans and therefore serve as a convenient benchmark. Their distinctive 




(41) 
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Function 


Image 


Template 


Temporal 
integration 


4>o = —~4>o + k • 9if(t, 6 2 ) 

T 


fa = --fa + k ■ Oi if(t, 9 2 ,i) 

T 


Adaptation 
to brightness 


No 


Ol,i = (00 - fahi + Aj,i 
7i 

Aj,i = —(fa -fa) 

T 


Adaptation 
to rotation 


No 


h.i = (A 2ji (t) + l)7r 

Aj,2 = 72Aj,3 ||0O - fa ||e 
Aj,3 = -72Ai,2||0O - fa\\e 


Detectors 
of similarity 


x = -axl + bxl + y - z + I 

+uo + <Po(t), 
y = c- dx 2 Q- y , 
z = e(s(x + x Q ) - z ), 


±i = -ax\ + bxf + yi - Zi + I 

+Ui + fa(t), 
yi = c - dxf - yi, 
Zi = e(s(xi + Xi) - Zi), 


Coupling 
function 


u = 7 (-(n + l)x + T,j^o x j) 


iH = 7 (-(«+ + Ej^^i) 



Table 4: Equations of the system for rotation and brightness-invariant template matching. 
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Figure 8: Template matching of Garner patterns of various complexity. The patterns (upper 
row) were rotated by tt/4 and had various intensity. Depending on the number of their 
rotational symmetries they induced different number of invariant sets in the system state 
space: two, four and eight respectively. The diagrams of corresponding phase plots are 
provided in the middle row. Estimates of the rotation angle as functions of time for different 
initial conditions are depicted in the third row. 

property is that their overall intensity does not vary from one pattern to another. At the 
same time their complexity in terms of the number of rotation and reflection symmetries can 
vary. In our example the first pattern has the highest complexity. The second has mid level of 
complexity (two symmetrical states) and the third has the lowest degree of complexity (four 
symmetrical states). The second row of Figure [8] illustrates the system dynamics involved in 
invariant template matching for these patterns. The diagrams represent phase plots of the 
successful node j (e.g. for the template subsystem in which the matching occurs). The third 
row contains trajectories of the estimates of the rotation angle 6j^{t). Each object induces 
various number of invariant sets in the template subsystems. The number of these invariant 
sets is inversely proportional to stimulus complexity. Hence, the higher the complexity 
the more time the system requires to converge to an attractor. Thus the time needed for 
recognition increases monotonically with the stimulus complexity. This is consistent with 
empirical results reported in many experimental studies, for instance [25] . 

An additional property of our system is that it is capable of reporting multiple represen- 
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tations of the same object. This is indicated by the dashed trajectories in Figure [SJ Even 
though the system parameters are chosen such that trajectories converge to an attractor, 
we can still observe meta-stable behavior. This is because the attractors in our system are 
of Milnor-type, which implies that trajectories starting in the vicinity of one attractor may 
actually belong to the basin of another attractor. Furthermore, it is even possible to tune 
the system in such a way that it will always switch from one representation to another. The 
latter property suggests that our simple system in Table H] can provide a simple model for 
visual perception, where spontaneous switching and perceptual multi-stability are commonly 
observed [2], [26] . 

5.5 Tracking disturbances in scanning microscopes 

We next consider the application of a the template-matching system with weakly attracting 
sets to a problem of realistic complexity. We applied our approach to the problem of tracking 
morphological changes in dendritic spines based on measurements received from a multipho- 
ton scanning microscope in vitro. A distinctive property of laser microscopy is that in order 
to "see" an object one needs, first, to inject it with a photo-sensitive dye (fluorophore). The 
particles of the fluorophore emit photons of light under external stimulation, thus illumi- 
nating an object from inside the tissue. Typical data from a two-photon microscope are 
provided in Figure [9P. 

We addressed the problem of how to register fast dynamical changes in spine geome- 
try after application of chemical stimulation. The measurements were performed on slices. 
Measurements of this kind suffer from effects of photobleaching and diffusion of the dye (see 
Figure [9]), and dependance of the scattering of the emitted light on the a-priori unknown po- 
sition of the object in the slice. On-line estimation and tracking the effects of photobleaching 
(intensity) and changes of the object position (blur) in the slice are therefore necessary. 

The measured signal is already a temporal sequence, which fits nicely to our approach. 
An inherent feature of scanning microscopy is that the object is measured using a sequence 
of scans along one- dimensional domains (see Figure El panel a). Hence the objects in this 
case are one-dimensional mappings, and the domain Q x is an interval Q x = [x m i n , x max ]. For 
the particular images we set x m i n = 1 and x max = 176, which corresponds to a scanning line 
of 176 pixels and 5.95 micro meters. In order to eliminate measurement noise we we consider 
the averaged data in the scanning line over n successive subsequent trials. 

The measured image, So, was chosen to be the averaged data along the scanning line over 



The images are provided by S. Grebenyuk, group of neuronal circuit mechanisms, RIKEN BSI 
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Figure 9: Typical images from the two-photon microscope. Panel a shows a dendrite; 
the domain of scanning (red line) is in the vicinity of two spines (small protrusions on the 
dendrite). Size of the domain is 5.95 micron, and speed of scanning, v s , is 1 pixel per 2 micro 
seconds. Panel b shows results of scanning as a function of time in the beginning (interval 
[Ti,T2]), in the middle of experiment (domain p^Tj), and in the end of the experiment 
(domain [T 5 ,T 6 ]). 

n successive subsequent trials. The template, Si, substituted the averaged measurements of 
the object at the initial time T\. 

Samples of data used to generate Si are provided in Figure [TU1 a. These correspond to 
the intensity of the emitted radiation from the object in the red part of the spectrum for 
the data shown in Figure [9j b, fragment 1. Measured objects, So, are the averaged samples 
of data at the time instants Tj ^ T\ (proportional to T s ). Focal distortions were simulated 
using conventional filters from Photoshop applied to Si. These fragments are provided in 
Figure [10J panels b and c. 

Because the sources of perturbation are the effects of photobleaching (affecting bright- 



a b c 




min 



Ti-T m T x T 2 -T m T 2 T 3 -T m T 3 

Figure 10: Data which has been used to generate the template, Si (panel a), and perturbed 
measurements So at time instants T2 and T3 (panels b and c respectively). 
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Figure 11: Trajectories e(t), 0i(t), ^(t) as functions of time. Black lines correspond to 
measurements in Fig. [TO], panel b. Blue lines correspond to the data in Fig. [TO] panel c. 
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Figure 12: Plots of the synchronization errors xo(t) — Xi(t) as a function of time. Panel a 
corresponds to the data depicted in Fig. [TO] b. Panel b corresponds to the measurements 
shown in Fig. [10] c. 



ness) and deviations in the object position in the slice (affecting scattering and leading to 
blurred images) the following model of uncertainty was used: 



9ifi(x,9 2 ,t)=9 1 [ e-^-^S^K, 
Jn x 

where x(t), the scanning trajectory in (1421) . is defined as: 



(42) 



x(t) 



t ^ X m ^- X Xr^ 



x(t (imax 3'min)); t ^> X 



max ^min 



, k s — 1. 



Figures \TT\ [T2] show the performance of our system <HM . fT2Tl) . fl3T?]) . fTTB]) in tracking 
focal/brightness perturbations for two measurements Sq. Figure [11] shows the tracking 
of unknown modelled perturbations in the images. Figure [12] shows the synchronization 
errors of the detection subsystem. Symbol t syQ denotes "synchronization time" spent in 
the vicinity of the invariant synchronization manifold. As follows from both figures, the 
system successfully tracks/reconstructs the estimates of unknown perturbations applied to 
the object (Fig. [TT]) . Coincidence detectors report synchrony only when the error between 
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the profiles of the template and object are is sufficiently small (Fig. [T2j) . Furthermore, the 
physical time required for recognition on the standard PC was less than 5 seconds. 

6 Conclusion 

We provided a principled solution to the problem of invariant template matching on the 
basis of temporal coding of spatial information. We considered the problem at the levels 
of mathematical analysis as well as implementation of specific recognition systems. Our 
analysis showed that a solution to the problem requires to abandon the traditional notion 
of attractor, in the Lyapunov sense, for defining a system's target set. As a substitute we 
proposed the concept of Milnor attracting sets. At the level of implementation we provided 
systems in which such attractors emerge as a result of external stimulation. These systems 
are endowed with mathematical rigor in the form of conditions sufficient for ensuring global 
convergence of trajectories to their target invariant sets. The results provided are normative 
in the sense that we require a minimal number of additional variables and consider as simple 
structures as possible. 

Even though the proposed system stems from theoretical considerations, it captures qual- 
itatively a wide range of phenomena observed in the literature on biological visual perception. 
These include multiple time scales for different modalities during adaptation [15], [Bj, [36] . 
switching and perceptual multi-stability j2], [26], perceptual ambiguity [43], [15], empirical 
observations in mental rotation [25] and decision time distributions [H] . This motivates our 
belief that present results may contribute to the further understanding of visual perception 
in biological systems, including humans. 

We demonstrated that the problem of invariant recognition can be solved by a simple 
system of ordinary differential equations with locally Lipschitz right-hand side. This result 
can be used as an existence proof for solving the problem of adaptive recognition by means 
of recurrent neural networks with fixed weights. Such systems are being used in various 
computational tasks [21] without any guarantee of a solution. We guarantee that solutions to 
realistic recognition problems (e.g. insuring invariance to rotation, blur, scaling, translation 
etc.) can be obtained with networks approximating our prototype system sufficiently well. 

Appendix 1 Optimality of sampled representations 

Consider an image S(x, y) and its quantized version S q obtained from S(x, y) by dividing 
domain fl x x Vl y into the union of finite number of subsets Vt x ^ x Vt y ^, Vt x = Uj x Q x> j, 
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Qy = uf y Q yi i. To each subset Q x j x £l y i a value is assigned, which can be thought of as the 
median value of S(x, y) over Q x j x fl yti . We represent S q as a function of indices S q (i,j) 
and assume that the value of S g (i,j) is quantized by a set of N s levels. 

Consider a system of sensors which are capable of measuring image S(x,y) instanta- 
neously over the given fc-union of subsets Q x j x Vt y ^. The system's cost can be naturally 



n order to measure the entire image at once 
so in the optimal case its cost C should 



defined in terms of its total number of sensors, 
the system must have at least N x N y /k sensors 3, 
equal C(k) = N x N y /k. 

We estimate the amount of information contained in this sampled representation of the 
image. The image is represented by an N x N y /k-tup\e of elements. Each element is assigned 
a value, say <Tj, from a set of N s levels with the given probability p(o~i). Hence the entropy 
of the representation is 

H ( k ) = -J2 p ^ log (jvT p ^)) = log (~T^) ~ Yl p ( a *) l °ZP( a ^ 

The entropy characterizes the informational content of a representation, and function 1/H(k) 
its ambiguity. 

Overall losses, Q(k), therefore can be defined as a weighted sum of costs, C(k), and 
ambiguity, 1/H(k): 

Q(k) = X 1 C(k) + X 2 l/H(k), Ai, A 2 e M >0 , k e [1, N x N y ] 

Function Q(k) is unimodal and increasing towards the boundaries of k: k = 1, k = N x N y . 
This implies that the minimum of Q(k) is achieved for some k = k* G (l,N x N y ). In other 
words, a representation is optimal only when it is sampled, e.g. induced by a finite, yet 
neither complete nor elementary, partition of the domain Q x x Vt y . 

Appendix 2 Proofs of the theorems 

Proof of Theorem^ We prove the theorem in three steps. First, we show that the solution of 
the extended system (T24|) . (j27j) . (!30|) is bounded. Second, we prove that there are constants 
p, b, e and time instant t' > such that the following holds for system solutions: 

WMt) ~ Mt)\U < e- p(t - to) \\Mto) - Mh)\\e + b\\9 2 - M^IUm V t > t > t' (43) 

Third, using this representation we invoke results from (our paper) and demonstrate that 
conclusions of the theorem follow. 



For simplicity wc assume that N x N y can be expressed as multiples of k. 
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1. Boundedness. To prove boundedness of solutions of the extended system in forward 
time let us first consider the difference ei(t) = <po(t) — <pi{t). According to ( 1241) . dynamics of 
ei{t) will be defined as 

ei = -\e, + k (Oj^t, 8 2 ) - OtMt, 0;, 2 )) + e(t) (44) 

Noticing that 

OiMt, e 2 ) - OiMt, e i>2 ) = faMt, e 2 ) - ej^t, e ij2 )\ + [OxUit, e i>2 ) - e lA f t (t, e ij2 )\ 

and denoting 5\ = 6*0 — 9\, S 2 (t, 81, 8 2 , 6^2) = 8ifi(t, 82) — 9ifi(t, 8 iy2 ) we can rewrite (jSj) as 
follows 

et = --et - S^kfiit, 8 i>2 )\ + 5 2 (t, 9 U 8 2 , 8 i>2 )k + e(t) (45) 

T 



Let us now write equations for Qi.\—Q\ in (1261) in differential form. To do so we differentiate 
variable 8 it \ — 9± = 8\ with respect to time, taking into account equations (14"4|) . ( l45l) : 

Si = -7i (S 1 [kf l (t,8 lt2 )}-5 2 (t,8 1 ,8 2 ,8 h2 )k-e(t)~) (46) 

Variable e(t) in (1461) is bounded according to ( l23l) . Let us show that S 2 (t,9i,9 2 ,8^ 2 ) is 
also bounded. First of all notice that the following positive definite function 

V x = 0.5 (A2 + A3) 

is not growing with time: 

V = A 272 A 3 ||0o(t) - (f>M\U ~ A372A2II0OW - 

Furthermore 



\ 2 (t) = r ■ sin (l2 

t 



^o(t) - (f)i(r)\\ e dT + ipo 



\ 3 (t) = r - cos 72 / ||^o(t) - <Mt)|Mt + y? , r,<p e 



to 



(47) 



Choosing initial conditions \ 2 (to) + Ag(to) = 1 ensures that r — 1. Hence, according to 
equation (1301) . variable 8 i;2 belongs to the interval [8 2tinin , 8 2>iaax ]. 
Consider variable 5 2 (t, 61, 8 2 , 8^ 2 ): 

5 2 (t, 6 X , 2 , 8 lt2 ) = 8 1 f i (t, 8 2 ) - Oiffa 8 h2 {t)) = 8, (ji(t, 8 2 ) - /<(*, 0, >2 (i)) (48) 
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Taking into account notational agreement and properties (J5J), ([H]), we conclude that the 
following estimate holds 

|* 2 (*, 01, fe, ^i, 2 )| < » 3 ) - ^i,2(*)| < ^xn^c^^l^ - ^<, 3 (*)| (49) 

Given that 9i, 2 (t) G [#2,min, #2,max] and using ( l49l) we can provide the following estimate for 
5 2 (t, 9 1 , 2 , 9 it2 ): 

\^2{t, 6l, 02, 9i,2)\ < ^l.max-D^l^.max — ^2,min| (50) 

Let us consider equality ( )46l) . According to condition 1) of the theorem, term 

a{t) = kfi{tA,2{t)) 
is nonnegative and bounded from below: 

a(t) = kfi(t, hfiit)) > kD 3 , V t > (51) 
Taking into account equations ( 1461) . ( I5TT) we can estimate \S±(t)\ as follows: 

Mt)\ < e" Ylf *o a(T)dT \5 1 (to)\+jie-' riS *o aiT)dT f e 7l/ 'o a(ri)dT1 |e(r) + 5 2 {r)k\dr (52) 

J t 

According to (B5j). (1511)) we have that for all* > t > 

\e{t) + 5 2 {t)k\ < ||e(r) + fc<y 2 (r)||oo,[to,t] < A + k0 ltmax DD 2 \6 2 ,max ^2,111111 

| = Mi (53) 

Furthermore 



f() 



7i V a W a(t )/ 7i-°3fc 



Taking into account fl52|) . (1531 . and fl5^|) we can obtain the following estimate: 

|5 1 (t)|<e-^- t «)|5 1 (t )| + ^ : (55) 

D 3 k 

Inequality (l55j) proves that <5i(£) is bounded. 

In order to complete this step of the proof it is sufficient to show that e^it) is bounded. 
This would automatically imply boundedness of 4>i(t), thus confirming boundedness of state 
of the extended system. To show boundedness of Ci(t) let us write the closed-form solution 
of (HU): 

rt 

(*-*o) * / T l 

ei(t) = e ~ei(t ) + e~~ e~ (5i(ri)a(ri) + kS 2 (ri) + e(ri)) dn (56) 



to 



Using (l53l) and (133]) we can derive that 



|e,(t)| < e-^| ei (to)| + M,r (l + ^ ) + (57) 
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where ei(t) is an exponentially decaying term: 

l _ e -(^-nkD 3 )(t-toY 



-7ifcD 3 (i-t ) 



|#i — Oi,l(t )\. 



(5* 



As follows from (|47|) . (I55j) . (1371) . (1581) . variables ej(t), ^,2(0 are bounded. Hence state 

of the extended system is bounded in forward time. 

2. Transformation. Let us now show that there exists a time instance if and constants 
p, c e ]R >0 such that the dynamics of e*(t) = <f>o(t) — (f>i(t) satisfies inequality (l43l) . In order 
to do so we first show that term 

in (H5l) can be estimated as 



\Si(t)kfi(t, e i>2 (t))\ < M 2 \9 2 - e ii2 (t)\ + A 2 + e 2 (t) 



(59) 



where M 2 , A 2 are positive constants and e 2 (t) is a function of time which converges to zero 
asymptotically with time. 

According to (1521) the following holds 

< e" 71 J *o a(T) " r |^i(^o)| + 7ie~ 71 J '° a(T)dT / e 71 a(n)(iri |e(r) + <5 2 (r)£;|dr 



to 



Taking into account (123|) . (15TT) we can conclude that 

|5i(t)| < e-T 1 * 03 ^-*") |*i(*o)l + + 7ie" 7l/t o° (T)dr / V 1 ''o Q{ri)dri |5 2 (r)£;|dr (60) 



to 



Substituting (H9j) into (1601) results in 



A 



fc£>. 



J t 



(61) 



Consider the following term in ( IBTj) : 



(62) 



to 



Integration of (15"2"|) by parts yields 



g 7l Jig a(n)o!ri 



#2 - ^(r)|dr = - I -^-e^^ a{T>aT \e 2 - 9 i>2 (t)\ 



to 



7i \ Oi(t) 



| #2 — ^,2(^o) | 
Qt(t ) 



71 J tn «( r ) 



g 71 a(n)dri | ^2 ~ ^2^)1 



(63) 



o7i /A a(r)tir. - 



likD; 



■e TlJ '» a ^|e 2 -e i)2 ({)| + 



e 7i Jt «(ri)<iTi 



3 Jf 



rf|02-^ 2 (r) 
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Given that 



= e 2Mn + g2 ' max 2 ^ 2 ' min (A 2 (t) + 1), 



we can estimate the derivative d|0 2 — 02,i(t)\/dt as follows: 



5*2 — ^i,2 (^) ^ #2,max — #2,min i, /,\ , mi 

< — o — 72-|w)-<M*)| ( 64 ) 



dt 2 

Notice that the value of \4>o(t) — <fii(t)\ = ei{t) in (164)) can be estimated according to (1571) as 

|0o(t) - <f>i{t)\ < M x t (l + ^) + »i{t), 
where //i(t) ~ ei(£) + ej(t)e ~ is an asymptotically decaying term. 



Hence, taking into account (jolj) . fl57|) . flSTl) . fl63|) . and fl64j) we may conclude that the 
following inequality holds 

ir ^ 6l,maxDD 2] i , s, A 72 9i max DD 2 9 2>m ax ~ ^2,min „ , A . . /,s 

< — ^— 1^ - M*)l + + — J5p 2 MlT ( V 1 + ^J + Mt) 

where /x(£) is asymptotically vanishing term. Therefore f!59p holds with the following values 
of M 2 and A 2 : 

_ k9 hmax DD 2 D 4 

1V1 2 — ~ 

6, m DD 2 D., /, .Wfc^-W An (65) 



A 2 = 2H 
7i 



AL> 4 



To finalize this step of the proof consider variable ei(t) for t G oo], £i > to- According 
to ([56]), gSD we have that 

|e*(t) | < | ei (ti) | + rM 2 ||0 2 - i)2 (t) ||oo, [tl , t ] + 

^A 2 (^1 - e + r||e 2 (t)||oo,[t 1 ,oo] (l - e ~J + (66) 

rke hmax DD 2 \\e 2 - M^IUim + tA (i - e -^) 

Regrouping terms in ( 1661) yields: 

| ei (t)| - r (A 2 + A + ||e 2 (t)|| o,[t 1 ,oo]) < e"^ 1 (|e<(*i)| - r (A 2 + A + ||e 2 (t) lU^oo])) 

+ r(Af 2 + k6 hmax DD 2 )\\e 2 - ^ 2 (t)|U, [tl)t] 

Denoting 

A / = r(A 2 + A+||e 2 (OIU,[ tl ,cx,]) (67) 

we can obtain 

\e t (t)\ - A' < e-^ 1 (| ei (tOI - A') + t(M 2 + ^ 1>maxJ DL> 2 )||0 2 - <|2 (*)lk[t llt ] 



(68) 

< e-^We^Wu + r(M 2 + k6 hmax DD 2 ) \\9 2 - i , 2 (t)|| oo , [tl , t] 
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Given that 

' h(t)|-A', |ej(t)| > A' 



0, |ei(t)| < A' 

and taking into account inequality (|68|) . we can conclude that 

\\ei(t)\\ A , < e-^Wetih)^ + r(M 2 + k9 ltmax DD 2 )\\e 2 - 9 i<2 {t) \U, [tlA (69) 

Because equations ( 1661) - ( 1691) hold for any t\ G (to, oo] and that 

limsup ||e 2 (t)||oo,[ti,oo] = 

for every 

e > r(A + A 2 ) 

there exists a time instant t' > to such that the following inequality is satisfied 



< e-^|| ei (ti)|| 6 + r(M 2 + ke 1>msx DD 2 )\\9 2 - M^lkk,*] (70) 



for all t > t\ > t'. This proves (T431 for p = -, b = t(M 2 + k9i tmax DD 2 ). Hence the second 
step of the proof is completed. 

3. Convergence. In order to prove convergence we employ the following result from [38] : 

Lemma 1 (Corollary 3 in [38J ) Consider the following interconnection of two systems: 

S a : \\x(t)\\ A < \\x(to)\\ A -Pt(t-t ) + c- ||Mt)IIoo.[M> x:M> ^M n 

/•* /•* (71) 

£„: / 7||x(r)||^dr</i(to)-M*)< / 7 l|x(r)||^dr, V t > t , t eR + 

J to J t 

where the systems S a , S w are forward- complete^, function (3 t : M>o — > IR>o strictly mono- 
tone and decreases to zero as t — > 00. Lei suppose that the following condition is satisfied 

j-c-G<l, (72) 

/or some d G (0, 1), k G (1, 00). 

Then there exists a set f2 7 of initial conditions corresponding to trajectories x(£), /i(t) 
sitc/i that 

limsup ||x(t)||^ < c-h(t ); h(t) G [0, /i(i )] V t > t 

t— >oo 



8 We say that a system is forward-complete iff its state is defined in forward time for all admissible 
inputs. For the system S a the inputs are functions h(t) from Loo\tQ,t]. For the system S w the inputs are 
locally-bounded in t functions x(i). 
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In particular, f2 7 contains the following domain 

h(t ). 



|x(to)IU< 



In order to apply Lemma [1] we need to further transform equations (|47j) . flTOjl and 

M*) = 2 , mm + g2 ' max ~ g2 ' min (A 2 (t) + 1) (73) 

into the form of equation flTTj) . First, we notice that for every 6*2 G [6*2,111111, #2,max] there always 
exists a real number A* G [—1,1] such that 

/) n 1 ^2, max — $2. min / * ,\ 
0*2 = 02,min H » ^ A 2 + 1 J 

Hence, denoting 

l a /T i 7 7~i 7~i \ ^2, max — ^2 min 

c = r(M 2 + k6 ljmax DD 2 ) 

and using ( 1701) we ascertain that the following holds for solutions of system ( 1241) . ( 1271) . ( T30T) : 

||ei(t)|| 6 < e-^lle^OU, + c||A* - A 2 (t) ||oo,[t 1)t ] (74) 

for e > r(A + A 2 ), and t>h> t'. 

Consider the difference A?; — \2{t). According to (1471) we have 

\X* 2 - X 2 (t)\ < \a* - [ 72||e,(r)|| £ - ^1, \* 2 = sin(0 (75) 
Jti 

Denoting 

h(t) = a* - 72||ei(r)|| e - ip (76) 
Jti 

and taking into account (1741) . we therefore obtain the following equations 

\\ei(t)\\ £ < e-^He^OH, + c||/i(*)||oo,[ti,t] 

ft (77) 
/i(ti) - h(t) = / 72 1| ei(r) || e dr 

Equations (177j) are a particular case of equations (ITT]) to which Lemma [1] applies. In system 
(177]) . however, function is defined as j3 t {t) = e~r . Hence 

/3 t " 1 (t) = -rln(t) 

Therefore, according to Lemma [JJ satisfying inequality 

h (s)]^T(( 1 + T^) +1 ) <1 (78) 
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for some k G (1, oo), d G (0, 1) ensures existence of initial conditions ej(ii), h(ti) such that 
h(t) is bounded. Given that 

min In (^-) ( [1 + — J + 1 J « 15.6886 < 16 

«;6(i,oo)„de(o,i) Vet/ fc — 1 \\ 1 — dJ J 

we can rewrite condition (17H1) in a more conservative, yet simpler form: 

Taking into account notations (|65i) . ( 1761) we can rewrite this inequality as follows: 

nn ft i /^2,max — #2,min X 

k9 1 ^ mSLK DD 2 1 + — 



\4ry - ' D 3 J \ 2 

Notice that because the function sin(-) is periodic, the value of a* in (175!) and, subse- 
quently the value of h(ti), can be chosen arbitrarily large. Hence for any finite ej(ti) and <^o 
there will always exist a* and such that variable h(t) is bounded. 

Taking into account that h(t) is monotone and bounded, we can conclude that according 
to the Bolzano- Weierstrass theorem function h(t) has a limit in [0, h(ti)]: 



3h* G [0, fifa)] : lim h(t) = h*. 



t— >oo 



This in turn implies that 



Therefore 



lim / 72 1 1 ei(r) || e dr = a* - ip - h* < oo 



-ft 

"2, max "2, min 



3 6*2 G [0 2 ,min, #2,max] = lim 6 i>2 (t) = #2,min H ! (sm((7* - <p - h*) + 1) = 6' 2 

t— >oo 2 

Moreover, because ||ej(£)|| e is uniformly continuous in t, convergence of ||ej(t)|| e to zero as 
t — > oo follows immediately from Barbalat's lemma. T7ie theorem is proven. 

Proof of Theorem^ The proof consists of three major steps. First, we show that single 
Hindmarsh-Rose oscillator is a semi-passive system with radially unbounded storage function 
[30] . In other words, system: 

x = —ax 3 + bx 2 + y — z + I + u 

y = c — dx 2 — y (79) 
z = e(s(x + Xo) — z), a, b, c, d, e, s > 
obeys the following inequality 

V(x(t), y(t), z(t)) - V(x(0), ?/(0), z(0)) < f x(r)u(r) - H(x(r), y{r), z{r))dr. (80) 

Jo 
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where function H(-) is non-negative outside a ball in M 3 , and function V is positive definite 
and radially unbounded. Second, similar to [30], we show that semi-passivity of ( 1791) implies 
that solutions of the coupled system (|T6l) are bounded. Third, for an arbitrary pair 
of the oscillators we present a nonnegative function such that properties fl35l) . fl36l) hold for 
sufficiently large values of 7. Then we use the comparison lemma [21] to complete the proof. 

1) Semi-passivity of the Hindmash-Rose oscillator. Let us consider the following class of 
functions V: 

V(x, y,z) = - {cix 2 + c 2 y 2 + c 3 z 2 ) 

Then showing existence of a function V from the above class which, in addition satisfies 
inequality 

V < xu - H(x,y,z), (81) 

where H is non-negative outside some ball in IR 3 , would imply semi-passivity of (|79|) . 
Consider the time-derivative of V: 

V(x, y, z) = —c\ax A — c 2 dx 2 y — c 2 y 2 + c\xy 

— c 3 ez 2 + (c 3 es — c\)xz + c\bx 3 + c\Ix + c 2 cy + c 3 esxqz + c\xu. (82) 

Let us rewrite (1821) such that the cross terms xy, xz and x 2 y are expressed in terms of the 
powers of x, y, z and their sums. In order to do this we employ the following three equalities: 

- c 2 y 2 + Cl xy = -c 2 X 2 y 2 - c 2 (l - A 2 ) (y - + ^ _ ^ x 2 (83) 

2 / x , 2 f -i \ \ ( c 3 es - ci \ 2 (c 3 es - Ci) 2 2 

-c s ez + (c 3 es-c 1 )xz = -c 3 eX 3 z - c 3 e(l - A3) [z - ^ _ ^ arj + ^ _ 

(84) 

4 7 2 4 ( 2 C 2^ \ 2 (c 2 C?) 2 2 

— Ciax — c 2 dx y = —CiaXiX — Ciall — Ai) x H ; r^V H ; r^V (85) 

y v ; V 2cia(l-Ai)V 4cia(l-Ai) y v ; 

In what follows we will assume that constants Ai, X 2 and A3 in fl83|) - fl85|) are chosen arbitrarily 
in the interval (0, 1): < Aj < 1, % — 1, 2, 3. 

Taking equalities (I83p - (j85p into account, we can rewrite the time derivative of V (equation 
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in the following form: 
V(x, y, z) = - Cl a(l - Ai) (x 2 + 2 fl _ A , yJ - c 2 (l - A 2 ) (y - _ x) + 



- C3£(1 - As) (* - 2 < y(l-A3) g ) 2 ~ C2 ( A2 " AcMl- + ° 2Cy+ 

c 3 eX 3 z 2 + c 3 esx z - ciaAix 4 + cibx i + (- — — — + 7^7- — -r^r)x 2 + cilx + cixw 

V4c 2 (l-A 2 ) 4c 3 e(l-A 3 )/ 



Our goal is to express the right-hand side of (I86I) in the following form: 

V <dxu + (M - H (x,y,z)), (87) 

where H (x, y, z) is a radially unbounded nonnegative function outside a ball in R 3 , and M 
is a constant. For this reason we select constants A 2 , c 2 in (l8Tj) as follows: 



Noticing that 

c 2 d 2 



c 2 (\ 



4cia(l - Ai" 



2/ 2 + c 2 q/ 



/, c 2 rf 2 > / 2ccia(l - Ax) \ 2 cic 2 c 2 a(l - Ai) , , 

° 2[ 2 A Cl a(l- X 1 ) ) V 4A 2Cl a(l- Xx)-c 2 d?) 4A 2Cl a(l - A x ) - c 2 d 2 1 ' 

- c 3 eX 3 z 2 + c 3 esx z = -c 3 eX 3 (z - + C3 ^ g ° (90) 

proves representation (!H7|) for any fixed x = const. In order to show that ([87)1 holds with 
respect to the complete set of variables, e.g. (x, y, z) we use the following sequence of 
equalities: 

4 1 ( c i 2 [c 3 ss — C1) 2 \ 2 

— CiaAiX + C\bx + ; — ■ H ; — ■ x + ciix = (see notations below) 

V4c 2 (l-A 2 ) 4c 3 £(l-A 3 )/ 

— a x 4 + aix 3 + a 2 x 2 + a 3 x + a 4 = 



— fro^ 4 — (x — &i) 4 + b 2 x 2 + 6 3 x + 64 
2 + G? )x 2 — <i (x — di 

&o(x 2 - e ) 2 - (x - bi) 4 - d (x - rfi)^ + ei (91) 



6 x 4 — (x — 6i) 4 + (6 2 + d )x 2 — d (x — di) 2 + d 2 



with 

ci 2 (c 3 es-ci) 2 

a = CiaAi, ai = c x b, a 2 = - — + — r 

4c 2 (l - A 2 ) 4c 3 e(l - A3) 

a 3 = C\l , 0,4 = bo = cio ~~ 1) ^1 = b 2 = a 2 + |cti 2 , 6 3 = a 3 — y^ai 3 , 64 = 04 + 2 ^gCti 4 

do = 1, di = 7rr> ^2 = &4 + di 2 d , e = &2 ~|~ ^° , e x = d 2 + 6 e 2 
2a 2o 

(92) 
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Notice that we want the value of bo in (1911) . (I92p be positive. Hence the value of 



should be greater than 1. This can be ensured by choosing the value of c\ in ( 18TT) to be 
sufficiently large. As a result of this choice, taking restrictions (1881) into account, we conclude 
that the value of c 2 in flHTl) must be sufficiently small, e.g. satisfy the following inequality: 

4aA 2 (l- Ai) 



c 2 < Ci- 



d 2 
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The value for do can be chosen arbitrarily, here do — 1. 
Time-derivative V can now be written as follows 

V(x, y, z) = - Cl a(l - Ai) (z 2 + ^ _ Xi) y) 

' C3£(1 - Aa >(* - 2^1 ~V 

( sxo\ 2 c 3 es 2 x 2 

- Cs£X < z -2xJ 

(\ c 2 tP x / 2ccio(l - Ai) \ 2 cic 2 c 2 a(l - Ai) 

C2 ^ 2 4c 1 a(l-A 1 )H y 4A 2Cl a(l- A x ) -c 2 d?) 4A 2Cl a(l - A x ) - c 2 d 2 

— frof^ 2 — e o) 2 — ( x ~ bi) 4 — do{x — di) 2 + ei + C\xu (93) 
It is straightforward to see that expression ( J93l) is of the form (JHTl) . where 



tf (z, y, z) = Cl a(l - Ai) (z 2 + ^ + ^(1 - A 3 ) (s - J^l - V 

+ C2(1 - A2) (" " 2c,fl Cl -A,l X ) 2 + " 



2 



2c 2 (l-A 2 ) / d V 2A 

, c 2 d 2 . ( 2cc 1 a{l-\i] 
+ c 2 (A 2 - — ) ly 



4cia(l-Ai) / \ 4A 2 cia(l - Ai) - c 2 d 2 
+ b (x 2 - e ) 2 + (x - bi) 4 + d (x - d]) 2 

M _ c 3 es 2 x 2 dc 2 c 2 a(l - Ai) g 
4A 3 4A 2 cia(l - Ai) - c 2 d 2 Cl 

Let us denote 

H 1 (x,y,z) = H -M 

and rewrite (jSTjl as 

V" < cixm — Hi(x, y, z) 
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Function Hi(x, y, z) is radially unbounded. Furthermore, it is non-negative outside a ball in 
IR 3 . Hence choosing 

V*(x,y,z) = —V(x,y,z) 

Cl 

we assure existence of (radially unbounded) positive definite V*(x, y, z) such that 

V* < xu - H ^ z \ (94) 
ci 

where Hi(x, y, z)jc\ is radially unbounded and non-negative outside a ball in R 3 . Thus, 
according to (1HTT) . semi-passivity of the Hindmarsh-Rose system follows. 

2) Boundedness of the solutions. We aim to prove that boundedness of 4>i(t), i G 
{0, . . . , n} implies boundedness of the state of the coupled system. Without loss of gen- 
erality we assume that 

||<Mt)||oo,[0,oo] < 

Let us denote 

Vi = V*{xi,yi,Zi), H Xi = —Hi(xi,yi,Zi). 

Cl 

Consider the following function 

Vs(x, y, Z) = p V i( X ii Vi' Z i)' ■ ( 95 ) 

where x = col(x , . . . , x n ), y = col(y , • • • , Vn), z = col(z , ...,z n ) and 

l r , j s-C, s>C 
P^ C ^ = { 0, s<C 

Function Vs is nonnegative for any CeR and, furthermore, is radially unbounded. Hence, 
its boundedness for some C G R implies boundedness of Xi, yi, Zi, i G {0, . . . , n}. 
Let us pick C G R such that interior of the domain 



Q c = {x,y,zGl | y~]Vi(xi,yi,Zi) < C} 



i=0 

contains the domain 

n 

J2H hi (xi,yi,Zi) - kxI < M u Mi G R >0 , « G R >0 

i=0 

where M, is an arbitrarily large and k is an arbitrary small positive constant. In other words 
the following implication holds: 

n n 

^ Vi(xi, y u z^ >C ^ ^ H lti (xi, y^ zi) - kx\ > Mj (96) 

i=0 i=0 
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Such C always exists because Hi ti (xi, y^ z^ — kx 2 can be expressed as a sum of a nonnegative 
quadratic form in Xi,yi,Zi and non-negative functions of the higher order plus a constant, 
and Vi(xi,yi, Zi) is a positive-definite quadratic form. 

Consider time-derivative of function Vs(x, y, z). According to (195]) . fl94|) it is zero for all 
x, y,z G Qc, and satisfies the following inequality otherwise: 

n n n 

Ve < ^XiUi - y ^Hi^(xi,y i ,z i ) = 7x r Ix + ^Xjfajt) - ^H^Xj, yj, Zj) 

i=0 i=0 i=0 

Using Gershgorin's circle theorem, we can conclude that 

n n 

A Kp^ii y%i ^i) 

i=0 i=0 

Rewriting 

X4i{t) = K Li - + KXI + ^4>Kt), « > 

leads to the following inequality 

n n , D 2 \ n f D 2 \ 

i=0 i=0 ^ ' i=0 ^ ' 

Hence, choosing the value of C such that Mi > D^/Ak in (|96|) we can ensure that 

U E <0 

This implies that Vs(x(t), y(t), z(t)) is not growing with time. Hence trajectories Xi(t), yi(t), 
Zi(t) in the coupled system are bounded. 

3) Convergence to a vicinity of the synchronization manifold. Consider the i-th and j-th 
oscillators in ffl6l) . z, j g {0, . . . , n}, i ^ j. Let us introduce the following function 

V = 0.5 (C x ( Xi - Xj f + C y ( yi - Vj f + C z ( Zi - z 3 ) 2 ) , (97) 

where C x , C y > are to be defined and C z = G x j (se). 
Its time-derivative can be expressed as follows: 

V = -C x { Xl - x,) 2 + ^ + | Xj)2 ~ Kxi + x,) + 7 (n + 1)) (98) 

+C x (yi - yj){xi - £,•) - C^x; - + x^)^ - yj) 

-C y (yi - yj) 2 - C z e(zi - zj) 2 + C x (xi - Xj)(<pi - 4>j) 

Consider the following term in (l98j) : 

C x (yi - Vj)(xi - Xj) - C y d(xi - Xj)(xi + Xj)(yi - y s ) - C y (yi - yj) 2 . 
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It can be written as follows: 

-(l-Ai-As)^-^) 2 , (99) 
where Ai, A 2 G M>o and A x + A 2 G (0, 1). Taking (I9"9~|) into account we rewrite (I9"8"j) as: 



7 < -CM - x 3 ) 2 ( a -^ + ^ + a(Xl | Xj) - + x 3 f (100) 



-6(xi + Xj) + y(n + 1) - J - - Zi+i? 

-C y (l -A x - A 2 )(y t - y f + C x {x { - Xj)(fa - fa) 



Let 



C„ 2aA 



Then 

\2 



y _ ZUU 2 

cZ = d 2 ' 



a / b\ 2 a f b\ 2 , JX d 2 



F < -Or(zi - xj) 2 - Si - - +- s r - + 7(n + 1) 



2 V a/ 2 V a/ v ' 8aAiA 2 a 
-(1 - Ax - A2)C w (j/i - ytf - C z e( Zi - z i+1 ) 2 + C x (x t - Xj){fa - fa). (101) 



Hence, choosing 

1 / d 2 



7 > = I — + b 2 

1 (n + l)aV8AiA 2 

we can ensure that the first term in fllOll) is non-positive. The minimal value of 7 ensuring 
this property can be calculated by minimizing the value 

1 

8A!A 2 

for all Ai, A 2 G M> : A x + A 2 < 1. This can be done by letting A 2 = r - Ai, r G (0, 1) 
and differentiating the term 1/ (8A x (r — Ai)) with respect to Ax- This leads to the following 
solution: A x = r/2, A 2 = r/2. Taking this into account we rewrite (jlOip as follows 

tV / ~, ( a ( b Y a ( b Y t ^ d2 & 

v < CA.n •o) 2 ^ -) +- 2 {*>--) + 

-(1 - r)C y ( yi - y-j) 2 - C z e{ Zi - z i+l ) 2 + C x (x t - Xj)(fa - fa). (102) 
Let 

7 = 7 ~~~T\ (4 + &2 ) + £ i» £ i e K >o- 
(n + l)a V 2 y 
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Alternatively, we can rewrite this as 

7 = ( * + b 2 ) + e 2 , re (0, 1), e 2 E R >0 

(n + l)a \2r J 

Hence, according to f 1 1 2 j) the following inequality holds: 

V < -C x e 2 {xi - xj) 2 - (1 - r)C y (yi - yj) 2 - C z e(zi - z i+1 ) 2 + C x (xi - Xj)(</>i - (f>j). 
Then denoting a = 2 minj^, e, (1 — r)} we obtain 

V < -aV + C x ( Xi - - <j) 3 ) (103) 

Consider the following differential equation 

v = -av + C x (xi - Xj)(<f)i - <f)j) (104) 
Its solution can be estimated as follows 

\v(t)\ < e- a ^\v(to)\ + e~ at f e aT C x ( Xl (r) - x,(r))(^(r) - <P,(r))dr 
for all t > t . Given that Xi(t), Xj(t) are bounded there exists a constant B such that 

\v(t)\ < e~^\v(t )\ + ^-\\Mr) ~ 0i(r)IU,[*o, t] 
Then applying comparison lemma (see, for example [21], page 102) we can conclude that 

V(t) < e-^V(to) + ^-\\Ur) - Hr)\\oo M . 
Hence, conclusion 2) of the theorem follows. The theorem is proven. 
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